Population Parameter Estimation

Estimation of Friberg-Karlsson model of drug-induced myelosuppression.




Specification of prior distributions


Setup R environment

rm(list = ls())
gc()

modelName <- "friberg3"
scriptName <- paste(modelName, "Rmd", sep = ".")
fitModel <- FALSE

## Relative paths assuming the working directory is the script directory
## containing this script
scriptDir <- getwd()
projectDir <- dirname(scriptDir)
dataDir <- file.path(projectDir, "data")
modelDir <- file.path(projectDir, "model")
outDir <- file.path(modelDir, modelName)
figDir <- file.path(projectDir, "deliv", "figure", modelName)
tabDir <- file.path(projectDir, "deliv", "table", modelName)
invisible(dir.create(figDir, recursive = TRUE))
## Warning in dir.create(figDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/figure/friberg3' already exists
invisible(dir.create(tabDir, recursive = TRUE))
## Warning in dir.create(tabDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/table/friberg3' already exists
.libPaths("lib")

library(metrumrg)
## Loading required package: reshape
## Loading required package: lattice
## Loading required package: XML
## Loading required package: MASS
## Loading required package: grid
## metrumrg 5.57
## enter "?metrumrg" for help
suppressMessages(library(rstan))
suppressMessages(library(bayesplot))
suppressMessages(library(tidyverse))
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(parallel)
library(knitr)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(qapply)
## Loading required package: rlecuyer
## 
## Attaching package: 'qapply'
## The following object is masked from 'package:metrumrg':
## 
##     qstat
library(PKPDmisc)
library(GGally)
## 
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
## 
##     nasa
## The following object is masked from 'package:metrumrg':
## 
##     wrap
knitr::opts_chunk$set(
  comment = '.', 
  fig.height = 5, 
  fig.width = 9,
  eval.after = 'fig.cap',
  message = FALSE,
  fig.path = file.path(figDir, paste(modelName, "_", sep = ""))
)

## Go back to default ggplot2 theme that was overridden by bayesplot
theme_set(theme_gray())

source(file.path(scriptDir, "tools", "stanTools2.R"))
source(file.path(scriptDir, "tools", "functions.R"))

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

RNGkind("L'Ecuyer-CMRG")
set.seed(11191951) 
set.seed(10271998) 

The NONMEM model stub

modelFile <- file.path(modelDir, paste(modelName, "stub.ctl", sep = ""))
modelText <- readLines(modelFile)
cat(modelText, sep = "\n")
. $SIZES LVR=50
. $PROBLEM RUN 1 neutropenia PKPD example with Friberg-Karlsson PK-PD
. 
. $INPUT C ID TIME EVID AMT CMT DVO DV WT DVID
. $DATA ../../../data/friberg.csv IGNORE=(C='C')
. 
. 
. $ABBR DECLARE INTEGER FIRST_WRITE_PAR
. $ABBR DECLARE INTEGER FIRST_WRITE_IPAR
. 
. $SUBROUTINES ADVAN13 TOL=6 OTHER = ../../extrasend.f90
. 
. $MODEL
. NCOMPARTMENTS=8
. COMP = (GUT, DEFDOSE)
. COMP = (CENT)
. COMP = (PERIPH)
. COMP = (PROL)
. COMP = (TRANS1)
. COMP = (TRANS2)
. COMP = (TRANS3)
. COMP = (CIRC)
. 
. $PK
. ; Request extra information for Bayesian analysis.
. ; An extra call will then be made for accepted samples 
. include '/opt/NONMEM/nm74/util/nonmem_reserved_general'
. BAYES_EXTRA_REQUEST=1
. 
. nThin = THETA(12)
. 
. VWT = LOG(WT/70) ; normalized to 70 kg adult
. 
. MU_1 = THETA(1) + 0.75*VWT       ; CL
. MU_2 = THETA(2) + 0.75*VWT       ; Q
. MU_3 = THETA(3) + VWT            ; V1
. MU_4 = THETA(4) + VWT            ; V2
. MU_5 = THETA(5)                  ; KA
. MU_6 = THETA(6)                  ; MTT
. MU_7 = THETA(7)                  ; CIRC0
. MU_8 = THETA(8)                  ; ALPHA 
. MU_9 = THETA(9)                  ; GAMMA
. MU_10 = THETA(10)                ; log(SIGMA) additive error PK
. MU_11 = THETA(11)                ; log(SIGMA) additive error PD
. 
. " CALL EXTRASEND()
. 
. CL  = EXP(MU_1 +ETA(1))
. Q   = EXP(MU_2 +ETA(2))
. V1  = EXP(MU_3 +ETA(3))
. V2  = EXP(MU_4 +ETA(4))
. KA  = EXP(MU_5 +ETA(5))
. MTT = EXP(MU_6 +ETA(6))
. CIRC0=EXP(MU_7 +ETA(7))
. ALPHA=EXP(MU_8 +ETA(8))
. 
. GAMMA= EXP(MU_9 + ETA(9))
. SIG1 = EXP(MU_10 + ETA(10)) ; Additive Error PK
. SIG2 = EXP(MU_11 + ETA(11)) ; Additive Error PD
. 
. ; intermediate calculations
. K10 = CL/V1
. K12 = Q/V1
. K21 = Q/V2
. KTR = 4/MTT
. 
. KPROL = KTR
. KCIRC = KTR
. 
. ; initial PK conditions
. A_0(1) = 0
. A_0(2) = 0
. A_0(3) = 0
. 
. ; initial PD conditions
. A_0(4) = CIRC0
. A_0(5) = CIRC0
. A_0(6) = CIRC0
. A_0(7) = CIRC0
. A_0(8) = CIRC0
. 
. 
. S2 = V1
. S3 = V2
. 
. 
. IF(BAYES_EXTRA==1 .AND. NIREC==1 .AND. NDREC==1 .AND. &
.      ITER_REPORT>0 .AND. &
.      ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN 
. IF(FIRST_WRITE_PAR==0) THEN
. " OPEN(unit=52,FILE='./par.txt') 
.     FIRST_WRITE_PAR=1
.   ENDIF
.   " WRITE(52,'(I12,1X,50(1X,1PG19.10E3))') ITER_REPORT, &
.   " THETA(1), THETA(2), THETA(3), THETA(4), THETA(5), THETA(6), THETA(7), &
.   " THETA(8), THETA(9), THETA(10), THETA(11), &  
.   " OMEGA(1,1), OMEGA(2,1), OMEGA(2,2), &
.   " OMEGA(3,1), OMEGA(3,2), OMEGA(3,3), &
.   " OMEGA(4,1), OMEGA(4,2), OMEGA(4,3), OMEGA(4,4), & 
.   " OMEGA(5,1), OMEGA(5,2), OMEGA(5,3), OMEGA(5,4), OMEGA(5,5), &
.   " OMEGA(6,1), OMEGA(6,2), OMEGA(6,3), OMEGA(6,4), OMEGA(6,5), OMEGA(6,6), & 
.   " OMEGA(7,1), OMEGA(7,2), OMEGA(7,3), OMEGA(7,4), OMEGA(7,5), OMEGA(7,6), OMEGA(7,7), & 
.   " OMEGA(8,1), OMEGA(8,2), OMEGA(8,3), OMEGA(8,4), OMEGA(8,5), OMEGA(8,6), OMEGA(8,7), OMEGA(8,8)
. ENDIF
. 
. IF(BAYES_EXTRA==1 .AND. NDREC==1 .AND. ITER_REPORT>0 .AND. &
.   ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN 
.   IF(FIRST_WRITE_IPAR==0) THEN
.     " OPEN(unit=50,FILE='./ipar'//TFI(PNM_NODE_NUMBER)//'.txt') 
.     FIRST_WRITE_IPAR=1
.   ENDIF
.   " WRITE(50,'(I12,1X,F14.0,13(1X,1PG12.5))') ITER_REPORT, ID, &
.   " ETA(1), ETA(2), ETA(3), ETA(4), ETA(5), ETA(6), ETA(7), ETA(8)
. ENDIF
. 
. 
. 
. $DES
. 
. DADT(1) = -KA*A(1) 
. DADT(2) =  KA*A(1) - (K10 + K12)*A(2) + K21*A(3)
. DADT(3) =  K12 *A(2) - K21*A(3) 
. 
. CONC = A(2)/V1 ; Drug concentration
. 
. EDRUG = ALPHA * CONC
. IF(EDRUG > 1.0) EDRUG = 1.0
. 
. DADT(4) = KPROL *A(4) *(1 - EDRUG) * ((CIRC0 / A(8)) **GAMMA) - KTR *A(4)
. DADT(5) = KTR * (A(4) - A(5))
. DADT(6) = KTR * (A(5) - A(6))
. DADT(7) = KTR * (A(6) - A(7))
. DADT(8) = KTR*A(7) - KCIRC*A(8)
. 
. 
. $ERROR (OBSERVATIONS ONLY)
. IPRED = LOG(F)
. IND = 0
. IF(DVID.EQ.2) IND = 1 
. YPK = IPRED+SIG1*EPS(1) ; Error for PK
. YPD = IPRED+SIG2*EPS(2) ; Error for PD
. Y = YPK*IND + YPD*(1-IND)
. 
. ; Code below into chains

Assemble initial estimates, priors and table specs

## create initial estimates
geninit <- function(){
    paste(c(
    "; Initial values of THETA",
    "$THETA",
    paste(rnorm(1, log(10),0.5), "     ; Log CL = THETA(1)"),
    paste(rnorm(1, log(15), 0.5), "     ; Log Q = THETA(2)"),
    paste(rnorm(1, log(35), 0.5), "   ; Log V1 = THETA(3)"),
    paste(rnorm(1, log(105), 0.5), "   ; Log V2= THETA(4) "),
    paste(rnorm(1, log(2), 0.5), "   ; Log Ka = THETA(5) "),
    paste(rnorm(1, log(125), 0.2), "   ; Log MTT = THETA(6) "),
    paste(rnorm(1, log(5), 0.2), "   ; Log CIRC0 = THETA(7) "),
    paste(rnorm(1, log(3E-4), 1), "   ; Log ALPHA = THETA(8) "),
    paste(rnorm(1, log(0.17), 0.2), "   ; Log GAMMA = THETA(9) "),
    paste(rnorm(1, log(0.2), 0.2), "   ; Log PK Add Error = THETA(10) "),
    paste(rnorm(1, log(0.2), 0.2), "   ; Log PD Add Error = THETA(11) "),
    "$OMEGA BLOCK(8) ;INITIAL values of OMEGAs",
    paste(exp(rnorm(1, log(0.25), 0.5))),
    paste(rep(rnorm(1, 0.001, 0.01),1)), paste(exp(rnorm(1, log(0.25), 0.5))),
    paste(rep(rnorm(1, 0.001, 0.01),2)), paste(exp(rnorm(1, log(0.25), 0.5))),
    paste(rep(rnorm(1, 0.001, 0.01),3)), paste(exp(rnorm(1, log(0.25), 0.5))),
    paste(rep(rnorm(1, 0.001, 0.01),4)), paste(exp(rnorm(1, log(0.25), 0.5))),
    paste(rep(rnorm(1, 0.001, 0.01),5)), paste(exp(rnorm(1, log(0.25), 0.5))),
    paste(rep(rnorm(1, 0.001, 0.01),6)), paste(exp(rnorm(1, log(0.25), 0.5))),
    paste(rep(rnorm(1, 0.001, 0.01),7)), paste(exp(rnorm(1, log(0.05), 0.5))),
    "$OMEGA DIAG(3) ;INITIAL values of unused OMEGAs",
    "0 FIX; GAMMA",
    "0 FIX; Add PK",
    "0 FIX; Add PD",
    "$SIGMA  ;Initial value of SIGMA",
    "1 FIX; Add PK",
    "1 FIX; Add PD"),
    sep = "\n")
}

## Set parameters of the prior distribution
priors <- paste(c("$PRIOR NWPRI",
                  "$THETAP          ; Prior information of THETAS",
                  "(2.3 FIX)     ; CL Log(10)  ",
                  "(2.7 FIX)     ; Q  Log(15)" ,
                  "(3.56 FIX)     ; V1 Log(35)",
                  "(4.7 FIX)    ; V2 Log(105)",
                  "(0.693 FIX)      ; KA log(2)",
                  "(4.83 FIX)    ; MTT Log(125)",
                  "(1.61 FIX)      ; CIRC0 Log(5) ",
                  "(-8.11 FIX)   ; ALPHA Log(3E-4)",
                  "(-1.78 FIX)   ; GAMMA Log(0.17)",
                  "(-1.60 FIX)   ; log Add error PK Log(0.20)",
                  "(-1.60 FIX)   ; log Add error PD Log(0.20)",
                  "$THETAPV BLOCK(11)     ;  variances for priors on THETAS (var-cov)",
                  "0.25   FIX          ; CL ",
                  "0.00   0.25            ; Q",
                  "0.00   0.00   0.25        ; V1",
                  "0.00   0.00   0.00   0.25      ; V2",
                  "0.00   0.00   0.00   0.00   0.25    ; KA ",
                  "0.00   0.00   0.00   0.00   0.00   0.04    ; MTT ",
                  "0.00   0.00   0.00   0.00   0.00   0.00   0.04    ; CIRC0 ",
                  "0.00   0.00   0.00   0.00   0.00   0.00   0.00   1   ; ALPHA ",
                  "0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00  0.04  ; GAMMA  ",
                  "0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00  0.00  1.00  ; Error PK  ",
                  "0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00  0.00  0.00  1.00  ; Error PD  ",
                  "$OMEGAP BLOCK(8) ;prior information for OMEGA",
                  "0.045    FIX          ; CL ",
                  "0.00   0.045            ; Q",
                  "0.00   0.00   0.045        ; V1",
                  "0.00   0.00   0.00   0.045      ; V2",
                  "0.00   0.00   0.00   0.00   0.045    ; KA",
                  "0.00   0.00   0.00   0.00   0.00   0.045    ; MTT",
                  "0.00   0.00   0.00   0.00   0.00   0.00   0.045    ; CIRC0",
                  "0.00   0.00   0.00   0.00   0.00   0.00   0.00  0.045   ; ALPHA",
                  "$OMEGAPD (11, FIXED)     ; df for OMEGA prior"
),
sep = "\n")

## Table statements
tables <- paste(c("$TABLE ID EVID TIME DV IPRED CWRES CWRESI NPDE WT",
                  "NOPRINT ONEHEADER FILE=./.tab",
                  "$TABLE ID WT CL Q V1 V2 KA MTT CIRC0 ALPHA ETAS(1:LAST)",
                  "NOPRINT ONEHEADER FILE=./par.tab"),
                sep = "\n")

Run NONMEM

nChains <- 4
nPost <- 250 ## Number of post-burn-in samples per chain after thinning
nBurn <- 250 ## Number of burn-in samples per chain after thinning
nThin <- 1


seed = sample(10000, size = nChains)

if(fitModel){
  if(!file.exists(file.path(modelDir, modelName))) 
    dir.create(file.path(modelDir, modelName))
  
  runChains <- lapply(1:nChains, runChain,
                      modelName = modelName, modelDir = modelDir, 
                      priors = priors, tables = tables,
                      nPost = nPost, nBurn = nBurn, nThin = nThin, 
                      seed = seed, print = 100,
                      OLKJDF = 0, OVARF = 1, AUTO = 2,
                      NUTS_DELTA = 0.8,
                      grid = TRUE, method = "NUTS", mode='not para')
  

    }

read in data file

dataFile <- "friberg.csv"
nmData <- read.csv(file.path(dataDir, dataFile), as.is = TRUE, na.strings = '.')

Get population parameters

chains <- 1:nChains
##chains <- c(1, 3, 4)

nChains2 <- length(chains)

## Read in posterior distributions by chain and add a column called chain
popParameters = map(chains, function(thisChain){
  param <- data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""), 
                                       "par.txt"), data.table = FALSE) 
  names(param) <- c("iteration", 
                    "THETA1", "THETA2", "THETA3", "THETA4", 
                    "THETA5", "THETA6", "THETA7", "THETA8",
                    "THETA9", "THETA10", "THETA11",
                    "OMEGA11", "OMEGA21", "OMEGA22","OMEGA13", "OMEGA23", "OMEGA33",
                    "OMEGA14", "OMEGA24", "OMEGA34","OMEGA44", "OMEGA15", "OMEGA25",
                    "OMEGA35", "OMEGA45", "OMEGA55","OMEGA16", "OMEGA26", "OMEGA36",
                    "OMEGA46", "OMEGA56", "OMEGA66","OMEGA17", "OMEGA27", "OMEGA37",
                    "OMEGA47", "OMEGA57", "OMEGA67","OMEGA77", "OMEGA18", "OMEGA28",
                    "OMEGA38", "OMEGA48", "OMEGA58","OMEGA68", "OMEGA78", "OMEGA88")
  
  param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
  mutate(sample = 1:n())

## Calculate more interpretable paraneters
popParameters2 <- popParameters %>%
  mutate(CLHat = exp(THETA1),
         QHat = exp(THETA2),
         V1Hat = exp(THETA3),
         V2Hat = exp(THETA4),
         KaHat = exp(THETA5),
         MTTHat = exp(THETA6),
         CIRC0Hat = exp(THETA7),
         ALPHAHat = exp(THETA8),
         GAMMAHat=exp(THETA9),
         SigmaPK = exp(THETA10),
         SigmaPD = exp(THETA11)) %>%
  select(starts_with("OMEGA"),
         CLHat, QHat, V1Hat, V2Hat, KaHat, 
         MTTHat, CIRC0Hat, ALPHAHat, GAMMAHat, SigmaPK,
         SigmaPD, chain, iteration, sample)

## Convert to 3-D array with dims = {iterations, chains, parameters}. 
popParArray <- array(double(nPost * nChains2 * (ncol(popParameters2) - 3)), 
                     dim = c(nPost, nChains2, ncol(popParameters2) - 3), 
                     dimnames =  list(NULL, NULL, setdiff(names(popParameters2), c("chain", "iteration", "sample"))))
for(iChain in 1:nChains2){
  popParArray[,iChain,] <- popParameters2 %>%
    filter(chain == chains[iChain],
           iteration > 0) %>%
    select(-chain, -iteration, -sample) %>%
    as.matrix
}

Get individual parameters

indParameters = map(chains, function(thisChain){
  files <- list.files(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = "")))
  iparFiles <- files[grepl("^ipar", files)]
  param <- map(iparFiles, function(thisFile){
    data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""), 
                                thisFile), data.table = FALSE)
  }) %>% 
    bind_rows %>% select(1:10)
  names(param) <- c("iteration", "ID", "ETA1", "ETA2", 
                    "ETA3","ETA4","ETA5","ETA6","ETA7","ETA8")
  param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
  arrange(chain, iteration, ID) %>%
  left_join(popParameters)

MCMC diagnostics and posterior distributions of parameters

options(bayesplot.base_size = 12,
        bayesplot.base_family = "sans")
color_scheme_set(scheme = "brightblue")
myTheme <- theme(text = element_text(size = 12), 
                 axis.text = element_text(size = 12))

plot_mcmcHistory <- mcmcHistory(popParArray,
                                 nParPerPage = 5, myTheme = myTheme)
plot_mcmcDensityByChain <- mcmcDensity(popParArray, 
                                        nParPerPage = 16, byChain = TRUE, 
                                        myTheme = theme(text = element_text(size = 12), 
                                                        axis.text = element_text(size = 10)))
plot_mcmcDensity <- mcmcDensity(popParArray, nParPerPage = 16, 
                                 myTheme = theme(text = element_text(size = 12), 
                                                 axis.text = element_text(size = 10)))

plot_mcmcHistory
. [[1]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

. 
. [[2]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

. 
. [[3]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

. 
. [[4]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

. 
. [[5]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

. 
. [[6]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

. 
. [[7]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

. 
. [[8]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

. 
. [[9]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

. 
. [[10]]
Friberg-Karlsson PK-PD example: MCMC trace plots.

Friberg-Karlsson PK-PD example: MCMC trace plots.

plot_mcmcDensityByChain
. [[1]]
Friberg-Karlsson PK-PD example: Univariate marginal densities by chain.

Friberg-Karlsson PK-PD example: Univariate marginal densities by chain.

. 
. [[2]]
Friberg-Karlsson PK-PD example: Univariate marginal densities by chain.

Friberg-Karlsson PK-PD example: Univariate marginal densities by chain.

. 
. [[3]]
Friberg-Karlsson PK-PD example: Univariate marginal densities by chain.

Friberg-Karlsson PK-PD example: Univariate marginal densities by chain.

plot_mcmcDensity
. [[1]]
Friberg-Karlsson PK-PD example: Univariate marginal densities.

Friberg-Karlsson PK-PD example: Univariate marginal densities.

. 
. [[2]]
Friberg-Karlsson PK-PD example: Univariate marginal densities.

Friberg-Karlsson PK-PD example: Univariate marginal densities.

. 
. [[3]]
Friberg-Karlsson PK-PD example: Univariate marginal densities.

Friberg-Karlsson PK-PD example: Univariate marginal densities.

#Not ideal when there are many parameters
#mcmc_pairs(popParArray[,,!grepl("^OMEGA", dimnames(popParArray)[[3]])])
#mcmc_pairs(popParArray[,,grepl("^OMEGA", dimnames(popParArray)[[3]])])

ggpairs(filter(popParameters2) %>% 
          select(CLHat:GAMMAHat), 
        diag = list(continuous = 'barDiag'), 
        upper = list(continuous = 'points'), 
        lower=list(continuous='points'))
Friberg-Karlsson PK-PD example: Correlation Plots.

Friberg-Karlsson PK-PD example: Correlation Plots.

ggpairs(filter(popParameters) %>% 
          select(OMEGA11,OMEGA22,OMEGA33,OMEGA44,
                 OMEGA55,OMEGA66,OMEGA77,OMEGA88), 
        diag = list(continuous = 'barDiag'), 
        upper = list(continuous = 'points'), 
        lower=list(continuous='points'))
Friberg-Karlsson PK-PD example: Correlation Plots.

Friberg-Karlsson PK-PD example: Correlation Plots.

caption <- paste("Friberg-Karlsson PK-PD example:",
      c(rep("MCMC trace plots.", length(plot_mcmcHistory)),
        rep("Univariate marginal densities by chain.",
            length(plot_mcmcDensityByChain)),
        rep("Univariate marginal densities.",
            length(plot_mcmcDensityByChain)),
        rep("Correlation Plots.", 2)))

ptable <- monitor(popParArray, 
                  warmup = 0, print = FALSE) %>% 
  as.matrix() %>%
  formatC(3) %>%
  as.data.frame
write.csv(ptable, file = file.path(tabDir, paste(modelName, 
                                                 "ParameterTable.csv", sep = "")))

ptable %>% 
  rename(SEmean = se_mean, SD = sd, pct2.5 = "2.5%", pct25 = "25%", median = "50%",
         pct75 = "75%", pct97.5 = "97.5%", Neff = "n_eff") %>%
  mutate(parameter = rownames(.), "95% CI" = paste("(", pct2.5, ", ", pct97.5, ")", 
                                                   sep = "")) %>%
  select(parameter, mean, SD, median, "95% CI", Neff, Rhat) %>%
  kable(caption = "Summary of model parameter estimates.") %>%
  kable_styling(bootstrap_options = "striped", full_width = F)
Summary of model parameter estimates.
parameter mean SD median 95% CI Neff Rhat
OMEGA11 0.135 0.0866 0.111 (0.0434, 0.371) 652 1
OMEGA21 -0.0336 0.0667 -0.027 (-0.191, 0.0861) 840 0.999
OMEGA22 0.169 0.104 0.142 (0.0546, 0.479) 597 1
OMEGA13 -0.00534 0.062 -0.00315 (-0.133, 0.106) 553 1
OMEGA23 -0.00751 0.0669 -0.00681 (-0.151, 0.13) 621 1
OMEGA33 0.11 0.0761 0.0866 (0.0359, 0.329) 713 1
OMEGA14 0.0525 0.0782 0.0384 (-0.0632, 0.25) 692 1
OMEGA24 -0.0724 0.09 -0.0541 (-0.31, 0.0596) 673 1
OMEGA34 -0.00232 0.0749 0.000888 (-0.168, 0.141) 725 0.999
OMEGA44 0.199 0.142 0.162 (0.0611, 0.585) 715 1
OMEGA15 -0.00232 0.0662 -0.00181 (-0.13, 0.12) 614 1
OMEGA25 0.00947 0.0646 0.00845 (-0.128, 0.153) 451 1
OMEGA35 0.0196 0.0603 0.0121 (-0.0775, 0.158) 753 1
OMEGA45 -0.0165 0.0776 -0.0132 (-0.18, 0.122) 532 1
OMEGA55 0.117 0.0833 0.0938 (0.0353, 0.32) 878 1
OMEGA16 -0.016 0.0526 -0.0139 (-0.12, 0.0849) 578 1
OMEGA26 0.00682 0.0558 0.00827 (-0.116, 0.124) 613 1
OMEGA36 0.00529 0.047 0.00345 (-0.0756, 0.105) 737 1
OMEGA46 -0.00318 0.0619 -0.00595 (-0.118, 0.135) 698 1.01
OMEGA56 -0.00348 0.0509 -0.00131 (-0.109, 0.0931) 891 1
OMEGA66 0.095 0.0668 0.0771 (0.0339, 0.277) 673 1
OMEGA17 0.0144 0.0574 0.0126 (-0.0901, 0.143) 684 1
OMEGA27 0.0354 0.0666 0.0255 (-0.0659, 0.199) 603 1.01
OMEGA37 -0.00556 0.0537 -0.00447 (-0.11, 0.101) 706 1
OMEGA47 -0.0211 0.0764 -0.0157 (-0.181, 0.1) 773 1.01
OMEGA57 0.0104 0.0552 0.00674 (-0.0949, 0.131) 624 1
OMEGA67 -0.014 0.0483 -0.011 (-0.108, 0.0711) 713 1
OMEGA77 0.121 0.0815 0.101 (0.0412, 0.339) 493 1.01
OMEGA18 0.0372 0.062 0.0287 (-0.0583, 0.167) 586 1
OMEGA28 -0.0206 0.0628 -0.0165 (-0.159, 0.104) 777 1
OMEGA38 0.00601 0.052 0.0029 (-0.0856, 0.124) 669 1
OMEGA48 0.0265 0.0672 0.0212 (-0.0977, 0.17) 662 1
OMEGA58 -0.00619 0.0583 -0.00455 (-0.132, 0.114) 759 1.01
OMEGA68 -0.0103 0.0488 -0.00807 (-0.0983, 0.0716) 427 1
OMEGA78 0.0148 0.0564 0.0117 (-0.0807, 0.134) 497 1
OMEGA88 0.122 0.0836 0.101 (0.0421, 0.34) 531 1.01
CLHat 11.2 1.57 11 (8.35, 14.6) 832 1
QHat 22.6 3.63 22.4 (15.9, 30.2) 792 1
V1Hat 29.2 5.51 29.1 (18.7, 40.5) 801 1
V2Hat 117 20.4 114 (82.9, 162) 833 0.999
KaHat 1.69 0.344 1.66 ( 1.1, 2.42) 712 1
MTTHat 118 13.1 118 (93.7, 145) 741 1
CIRC0Hat 5.75 0.686 5.73 (4.46, 7.15) 956 1
ALPHAHat 0.000323 4.72e-05 0.000321 (0.000236, 0.000426) 780 1
GAMMAHat 0.151 0.013 0.15 (0.127, 0.177) 1.05e+03 1.01
SigmaPK 0.1 0.00442 0.1 (0.092, 0.109) 1.06e+03 1
SigmaPD 0.106 0.00678 0.106 (0.0941, 0.121) 908 1

Construct design data frame

## For now just use the NONMEM data set

design <- nmData %>%
  mutate(ID = id, WT = weight) %>%
  filter(is.na(C))

mrgsolve model code for posterior predictions

library(mrgsolve)

indCode =
  '
$PARAM 
THETA1 = 2.3
THETA2 = 2.7
THETA3 = 3.56
THETA4 = 4.7
THETA5 = 0.693
THETA6 = 4.43
THETA7 = 2.33
THETA8 = -1.78
THETA9 = -8.11
THETA10 = -1
THETA11 = -1
WT = 70
ETA1 = 0
ETA2 = 0
ETA3 = 0
ETA4 = 0
ETA5 = 0
ETA6 = 0
ETA7 = 0
ETA8 = 0

$CMT GUT, CENT, PERIPH, 
PROL, TRANS1, TRANS2, TRANS3, CIRC

$MAIN

double VWT = log(WT/70); 
double CL = exp(THETA1 + (VWT * 0.75) + ETA1 + ETA(1)); 
double Q = exp(THETA2 + (VWT * 0.75) + ETA2 + ETA(2));
double V1 = exp(THETA3 + VWT + ETA3 + ETA(3));
double V2 = exp(THETA4 + VWT + ETA4 + ETA(4));
double KA = exp(THETA5 + ETA5 + ETA(5));
double MTT = exp(THETA6 + ETA6 + ETA(6));
double CIRC0 = exp(THETA7 + ETA7 + ETA(7));
double ALPHA = exp(THETA8 + ETA8 + ETA(8));
double GAMMA = exp(THETA9);
double SIG1 = exp(THETA10);
double SIG2 = exp(THETA11);

// intermediate calculations
double K10 = CL/V1;
double K12 = Q/V1;
double K21 = Q/V2;
double KTR = 4/MTT;
double KPROL = KTR;
double KCIRC = KTR;

PROL_0 = CIRC0;
TRANS1_0 = CIRC0;
TRANS2_0 = CIRC0;
TRANS3_0 = CIRC0;
CIRC_0 = CIRC0;

$ODE
dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - K10*CENT - K12*CENT + K21*PERIPH;
dxdt_PERIPH = K12*CENT - K21*PERIPH;

double CONC = std::max(0.0,CENT/V1);
double EDRUG = CONC*ALPHA;
if(EDRUG > 1) EDRUG = 1;

dxdt_PROL = KPROL* PROL *(1 - EDRUG) * (pow((CIRC0 / CIRC), GAMMA)) - KTR *PROL;
dxdt_TRANS1 = KTR * (PROL - TRANS1);
dxdt_TRANS2 = KTR * (TRANS1 - TRANS2);
dxdt_TRANS3 = KTR * (TRANS2 - TRANS3);
dxdt_CIRC =  KTR*TRANS3 - KCIRC*CIRC;

$SIGMA 1 1

$OMEGA @block
0 
0 0
0 0 0
0 0 0 0
0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0 0

$TABLE
double cPredHat = std::max(0.000001,CENT/V1);
double lcPred = log(cPredHat) + SIG1*EPS(1);
double cPred = exp(lcPred);
double nPredHat = CIRC;
double lnPred = log(CIRC) + SIG2*EPS(2);
double nPred = exp(lnPred);

$CAPTURE lcPred lnPred cPred nPred
'

## Compile model
friMod <- mcode("FRImodel", indCode)

if(FALSE){
  ## Test runs
  out <- friMod %>%
    data_set(design) %>%
    carry.out(DV) %>%
    mrgsim %>% as.data.frame()
  
  ggplot(data=out, aes(x=time, y=cPred)) + geom_point(aes(group=ID, color=as.factor(ID)))
  
}

Individual predictions

samples <- unique(indParameters$sample)
set.seed(10301987)
indSim1 <- mclapply(samples,
                    function(thisSample){
                      friMod %>%
                        data_set(design) %>%
                        idata_set(indParameters %>% 
                                filter(sample == thisSample)) %>%
                        carry_out(sample, ID, time, DV, evid, DVID) %>%
                        mrgsim() %>%
                        as.tbl
                    }) %>%
  bind_rows

Population predictions

samples <- unique(popParameters$sample)
set.seed(10301987)
popSim1 <- mclapply(samples,
                    function(thisSample){
                      friMod %>%
                        data_set(design) %>%
                        param(popParameters %>% 
                                filter(sample == thisSample)) %>%
                        omat(omat(as_bmat(popParameters %>% 
                                filter(sample == thisSample),"OMEGA"))) %>%
                        carry_out(sample, ID, time, DV, evid, DVID) %>%
                        mrgsim %>%
                        as.tbl
                    }) %>%
  bind_rows

Posterior predictive distributions

## Combined plots

indPredC <- indSim1 %>% filter(evid==0, DVID==2) %>%
  group_by(ID, time, evid) %>%
  summarize(lbInd = quantile(cPred, probs = 0.05, na.rm = TRUE),
            medianInd = quantile(cPred, probs = 0.5, na.rm = TRUE),
            ubInd = quantile(cPred, probs = 0.95, na.rm = TRUE))

indPredN <- indSim1 %>% filter(evid==0, DVID==1) %>%
  group_by(ID, time, evid) %>%
  summarize(lbInd = quantile(nPred, probs = 0.05, na.rm = TRUE),
            medianInd = quantile(nPred, probs = 0.5, na.rm = TRUE),
            ubInd = quantile(nPred, probs = 0.95, na.rm = TRUE))

popPredC <- popSim1 %>% filter(evid==0, DVID==2) %>%
  group_by(ID, time, evid) %>%
  summarize(lbPop = quantile(cPred, probs = 0.05, na.rm = TRUE),
            medianPop = quantile(cPred, probs = 0.5, na.rm = TRUE),
            ubPop = quantile(cPred, probs = 0.95, na.rm = TRUE))

popPredN <- popSim1 %>% filter(evid==0, DVID==1) %>%
  group_by(ID, time, evid) %>%
  summarize(lbPop = quantile(nPred, probs = 0.05, na.rm = TRUE),
            medianPop = quantile(nPred, probs = 0.5, na.rm = TRUE),
            ubPop = quantile(nPred, probs = 0.95, na.rm = TRUE))

predAllC <- indPredC %>%
  left_join(popPredC) %>%
  left_join(design %>% filter(evid == 0, DVID==2)) %>%
  mutate(DV = suppressWarnings(as.numeric(DV)))

predAllN <- indPredN %>%
  left_join(popPredN) %>%
  left_join(design %>% filter(evid == 0, DVID==1)) %>%
  mutate(DV = suppressWarnings(as.numeric(DV)))

  p1 <- ggplot(predAllC, aes(x = time, y = DV))
  p2 <- p1 + 
    geom_line(aes(x = time, y = medianPop, 
                 color = "population")) +
    geom_ribbon(aes(ymin = lbPop, ymax = ubPop, 
                   fill = "population"), 
               alpha = 0.25) +
    geom_line(aes(x = time, y = medianInd, 
                  color = "individual")) +
    geom_ribbon(aes(ymin = lbInd, ymax = ubInd, 
                    fill = "individual"), 
                alpha = 0.25) +
    scale_color_brewer(name  ="prediction",
                       breaks=c("individual", "population"),
                       palette = "Set1") +
    scale_fill_brewer(name  ="prediction",
                      breaks=c("individual", "population"),
                      palette = "Set1")
  p3 <- p2 + geom_point() +
    labs(title = 'PK Model',
         x = "Time (h)",
         y = "PK Conc (ng/mL)") +
    theme(text = element_text(size = 12),
          axis.text = element_text(size = 12),
##          legend.position = c(0.8, 0.25),
          strip.text = element_text(size = 8)) +
    facet_wrap(~ ID, scales = "free_y")
  
  ####PD Plots#####

  p1b <- ggplot(predAllN, aes(x = time, y = DV))
  p2b <- p1b + 
    geom_line(aes(x = time, y = medianPop, 
                 color = "population")) +
    geom_ribbon(aes(ymin = lbPop, ymax = ubPop, 
                   fill = "population"), 
               alpha = 0.25) +
    geom_line(aes(x = time, y = medianInd, 
                  color = "individual")) +
    geom_ribbon(aes(ymin = lbInd, ymax = ubInd, 
                    fill = "individual"), 
                alpha = 0.25) +
    scale_color_brewer(name  ="prediction",
                       breaks=c("individual", "population"),
                       palette = "Set1") +
    scale_fill_brewer(name  ="prediction",
                      breaks=c("individual", "population"),
                      palette = "Set1")
  p3b <- p2b + geom_point() +
    labs(title = 'PD Model',
         x = "Time (h)",
         y = "ANC (10^9/L)") +
    theme(text = element_text(size = 12),
          axis.text = element_text(size = 12),
##          legend.position = c(0.8, 0.25),
          strip.text = element_text(size = 8)) +
    facet_wrap(~ ID, scales = "free_y")
  
  ####
  
  p3

  p3b

sessionInfo()
. R version 3.5.1 (2018-07-02)
. Platform: x86_64-pc-linux-gnu (64-bit)
. Running under: Ubuntu 14.04.5 LTS
. 
. Matrix products: default
. BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
. LAPACK: /usr/lib/lapack/liblapack.so.3.0
. 
. locale:
.  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
.  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
.  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
.  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
.  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
. [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
. 
. attached base packages:
. [1] parallel  grid      stats     graphics  grDevices utils     datasets 
. [8] methods   base     
. 
. other attached packages:
.  [1] mrgsolve_0.9.2        GGally_1.4.0          PKPDmisc_2.1.1       
.  [4] qapply_1.39           rlecuyer_0.3-4        kableExtra_1.1.0     
.  [7] knitr_1.23            gridExtra_2.3         forcats_0.4.0        
. [10] stringr_1.4.0         dplyr_0.8.3           purrr_0.3.2          
. [13] readr_1.3.1           tidyr_0.8.3           tibble_2.1.3         
. [16] tidyverse_1.2.1       bayesplot_1.7.0       rstan_2.19.2         
. [19] ggplot2_3.2.0         StanHeaders_2.18.1-10 metrumrg_5.57        
. [22] MASS_7.3-51.1         XML_3.98-1.20         lattice_0.20-38      
. [25] reshape_0.8.8        
. 
. loaded via a namespace (and not attached):
.  [1] httr_1.4.0                jsonlite_1.6             
.  [3] viridisLite_0.3.0         modelr_0.1.4             
.  [5] assertthat_0.2.1          highr_0.8                
.  [7] stats4_3.5.1              cellranger_1.1.0         
.  [9] yaml_2.2.0                pillar_1.4.2             
. [11] backports_1.1.4           glue_1.3.1               
. [13] digest_0.6.20             RColorBrewer_1.1-2       
. [15] rvest_0.3.4               colorspace_1.4-1         
. [17] htmltools_0.3.6           plyr_1.8.4               
. [19] pkgconfig_2.0.2           broom_0.5.2              
. [21] haven_2.1.1               scales_1.0.0             
. [23] webshot_0.5.1             processx_3.4.1           
. [25] generics_0.0.2            withr_2.1.2              
. [27] lazyeval_0.2.2            cli_1.1.0                
. [29] magrittr_1.5              crayon_1.3.4             
. [31] readxl_1.3.1              evaluate_0.14            
. [33] ps_1.3.0                  nlme_3.1-137             
. [35] RcppArmadillo_0.9.600.4.0 xml2_1.2.1               
. [37] pkgbuild_1.0.3            data.table_1.12.2        
. [39] tools_3.5.1               loo_2.1.0                
. [41] prettyunits_1.0.2         hms_0.5.0                
. [43] matrixStats_0.54.0        munsell_0.5.0            
. [45] callr_3.3.1               compiler_3.5.1           
. [47] rlang_0.4.0               ggridges_0.5.1           
. [49] rstudioapi_0.10           labeling_0.3             
. [51] rmarkdown_1.14            gtable_0.3.0             
. [53] inline_0.3.15             reshape2_1.4.3           
. [55] R6_2.4.0                  lubridate_1.7.4          
. [57] zeallot_0.1.0             stringi_1.4.3            
. [59] Rcpp_1.0.2                vctrs_0.2.0              
. [61] tidyselect_0.2.5          xfun_0.8